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Abstract 


A program  is  described  which  computes  Schwari-Christoh  'l  transformations 
that  map  the  unit  disk  conformally  onto  the  interior  of  a bounded  or  unbounded 
polygon  in  the  complex  plane.  The  inverse  map  is  also  computed.  The  computa- 
tional  problem  is  approached  by  setting  up  a nonlinear  system  of  equations  whose 
unknowns  are  essentially  the  "accessory  parameters"  ifc.  This  system  is  then  solved 
with  a packaged  subroutine. 

New  features  of  this  work  include  the  evaluation  of  integrals  within  the  disk 
rather  than  along  the  boundary,  making  possible  the  treatment  of  unbounded 
polygons;  the  use  of  a compound  form  of  Gauss-Jacobi  quadrature  to  evaluate  the 
Schwarz-ChristofTcl  integral,  making  possible  high  accuracy  at  reasonable  cost; 
and  the  elimination  of  constraints  in  the  nonlinear  system  by  a simple  change  of 
variables. 

Schwarz-Christoffel  transformations  may  be  applied  to  solve  the  Laplace  and 
Poisson  equations  and  related  problems  in  two-dimensional  domains  with  irregular 
or  unbounded  (but  not  curved  or  multiply  connected)  geometries.  Computational 
examples  are  presented.  The  time  required  to  solve  the  mapping  problem  is  roughly 
proportional  to  N3,  where  N is  the  number  of  vertices  of  the  polygon.  A typical 
set  of  computations  to  8-place  accuracy  with  N < 10  takes  1 to  10  seconds  on 
an  IBM  370/168. 
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I.  INTRODUCTION 


1.  Conformal  mapping  and  its  applications 

One  of  the  classical  applications  of  complex  analysis  is  conformal  map* 
ping:  the  mapping  of  one  open  region  in  the  complex  plane  C onto  another 
by  a function  which  is  analytic  and  one-to-one  and  has  a nonzero  deriva- 
tive everywhere.  Such  a map  preserves  angles  between  intersecting  arcs  in 
the  domain  and  image  regions;  hence  the  name  conformal.  The  Riemann 
Mapping  Theorem  asserts  that  any  simply  connected  region  in  the  plane 
which  is  not  all  of  C can  be  mapped  in  this  way  onto  any  other  such 
region.  The  theorem  does  not  say  what  this  mapping  may  look  like,  however, 
and  the  determination  of  particular  conformal  maps  for  particular  mapping 
problems  has  been  an  active  problem  since  at  least  1850. 


The  usefulness  of  conformal  mapping  for  applied  problems  stems  from 
the  fact  that  the  Laplacian  operator  transforms  in  a simple  way  under  a 
conformal  map.  Let  /: O -+C  map  a region  0,  in  the  s-plane  conformally 
onto  a region  Qw  in  the  u>-plane,  and  let  A,  and  Au,  denote  the  Laplacian 
operators  ab  + iy*  and  £?  + respectively,  where  s —»  x -+-  iy  and 
u>  = u -f-  »'v.  Then  we  may  easily  show, 

ArfM  - (i.i) 

for  ^:fl,  — »R  suitably  differentiable.  A conformal  map  has  |/f(*)|  > 0 every- 
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where;  thus  from  (1.1)  it  follows  that  if  <fr[z)  is  the  solution  to  the  Laplace 
equation  A,^  = 0 in  ft,,  subject  to  Dirichlet  boundary  conditions  <f>(z)  — 
g{z)  on  the  boundary  T, , then  ip(w)  = *(u>))  is  a solution  to  the  Laplace 

equation  A w\ft  = 0 in  the  image  region  ft*  = /(ft*),  subject  to  the  image 
boundary  conditions  ^(w)  = g(f~ l[w))  on  the  boundary  = /(T,).  (We 
have  assumed  that  / maps  T,  bijectively  onto  the  boundary  of  ftu,.  This  is 
not  always  true,  but  it  is  true  if  both  regions  are  bounded  by  Jordan  curves. 
See  [Henrici,  1974],  Thm.  5.10e.) 

More  generally,  from  (1.1)  we  can  see  that  Poisson's  equation,  A /f>(z)  = 
p[z ),  transforms  under  a conformal  transformation  into  a Poisson  equation 
in  the  id-plane  with  altered  right  hand  6ide: 

a„^w)  = i/'(r‘H)r1Kr‘N)-  (i.») 

Furthermore,  more  general  boundary  conditions  than  Dirichlet  also  trans- 
form in  a simple  way.  For  example,  the  Neumann  condition  — M*)> 

where  jf-  is  a normal  derivative  in  the  2-plane,  transforms  to  = 

{ w )).  We  do  not  pursue  6uch  possibilities  further  here; 
for  a systematic  treatment  see  chapter  VI  of  [Kantorovich  & Krylov,  1958]. 
Some  computed  examples  are  given  in  Section  V. 

Traditionally,  conformal  mapping  has  been  applied  mo6t  often  in  two 
areas.  One  is  plane  electrostatics,  where  the  electrostatic  potential  <p  satisfies 
Laplace's  equation.  The  other  i6  irrotational,  nonviscous  fluid  flow  in  the 
plane,  which  may  be  described  in  terms  of  a velocity  potential  <p  that  also 
satisfies  Laplace's  equation. 


2.  The  Schwars-ChristofTel  transformation 

The  problem  of  mapping  one  complex  region  conformally  onto  another 
is  in  general  very  difficult,  but  for  the  special  case  of  polygonal  regions  it 
can  be  greatly  simplified.  Suppose  that  we  seek  a conformal  map  from  the 
unit  disk  in  the  s-plane  to  the  interior  of  a polygon  P in  the  tu-plane  whose 
vertices  are  wi,...,  numbered  in  counterclockwise  order.  For  each  k, 
denote  by  the  exterior  angle  of  P at  u >*: 


For  any  polygon  we  have  a simple  relationship  among  the  numbers  /?*: 


(1.3) 

If  to*  is  a finite  vertex,  we  have  — 1 < /?*  < 1.  We  need  not  require,  however, 
that  P be  bounded.  It  may  have  a number  of  vertices  at  complex  infinity,  and 
the  exterior  angles  corresponding  to  these  may  fall  anywhere  in  the  range 
1 < Pk  < 3.  Such  angles  are  defined  to  be  equal  to  2 x minus  the  external 
angle  formed  in  the  plane  by  the  intersection  of  the  two  sides  involved,  if 
they  are  extended  back  away  from  infinity.  The  following  example  should 
illustrate  what  is  meant  by  various  values  of  /?*:  it  is  a polygon  with  five 
vertices  to*  (in  this  case  toi  = to*),  with  corresponding  values  (A,...,  ft)  = 

(i»  1): 


As  always,  (1.3)  holds  for  this  example. 

Let  us  now  pick  at  random  N points  afc  (“prevertices")  in  counterclock- 
wise order  around  the  unit  circle  and  two  complex  constants  C and  toc,  and 
consider  the  Schwars-ChristofTel  formula: 


The  quantities  (1  — r'Azfc)  always  lie  in  the  dUk  |u> — 1|  < 1 for  \z\  < 1. 
Therefore,  if  we  choose  a branch  of  log(z)  with  a branch  cut  on  the  negative 
real  axis  by  means  of  which  to  define  the  powers  in  (1.4),  w(z)  defines  an 
analytic  function  of  z in  the  disk  \z\  < 1,  continuous  on  \z\  < 1 except 
possibly  at  the  vertices  z*. 

The  Schwarz-ChristofTel  formula  is  chosen  so  as  to  force  the  image  of 
the  unit  disk  to  have  corners  in  it  with  the  desired  exterior  angles  /?**.  It 
i6  not  hard  to  see  from  (1.4)  that  at  each  point  z*,  the  image  u>(z)  must 
turn  a corner  of  precisely  this  angle.  This  is  in  keeping  with  our  purpose  of 
mapping  the  disk  onto  the  interior  of  P.  What  the  map  will  in  general  fail 
to  do  is  to  reproduce  the  lengths  of  sides  of  P correctly,  and  to  be  a one- 
to-one  correspondence.  For  a suitable  choice  of  parameters  {z/t},  C,  and  u»c, 
the  image  under  / of  the  unit  disk  might  be,  for  example, 


Only  the  angles  are  guaranteed  to  come  out  right. 

The  variables  zi,  ...,zn,  C,  and  uv  are  the  accessory  parameters  of  the 
Schwarz-ChristofTel  mapping  problem.  Our  first  problem — the  parameter 
problem — is  to  determine  values  of  the  accessory  parameters  so  that  the 
lengths  of  sides  of  the  image  polygon  do  come  out  right.  The  central  theorem 
of  Schwarz-ChristofTel  transformations  asserts  that  there  always  exists  such 
a set  of  accessory  parameters: 


Theorem  1 (Schwars-Christoffel  transformation).  Let  D be  a simply 
connected  region  in  the  complex  plane  bounded  by  a polygon  P with  vertices 
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z\,...,zh  and  exterior  angles  xflk,  where  — 1 < 1 if  Zk  is  finite  and 

1 < ft  < 3 if  Zk  — oo.  Then  there  exists  an  analytic  function  mapping  the 
unit  disk  in  the  complex  plane  conformally  ontoD,  and  every  such  function 
may  be  written  in  the  form  (1.4). 

Proofs  [Henri.ci,  1974],  Thm.  5.12e. 


In  fact,  for  any  given  polygon  there  is  not  just  one  but  infinitely  many 
such  conformal  mappings.  To  determine  the  map  uniquely  we  may  fix  ex- 
actly three  points  Zk  at  will,  or  fix  one  point  Zk  and  also  fix  the  complex 
value  wc,  or  (as  in  a standard  proof  of  the  Riemann  mapping  theorem)  fix 
wc  and  the  argument  of  the  derivative  /'( 0). 

The  simplicity  of  the  explicit  formula  (1.4)  is  attractive.  But  because 
the  problem  of  determining  the  accessory  parameters  is  intractable  analyti- 
cally, applications  of  it  have  almost  always  been  restricted  to  problems 
simplified  by  having  very  few  vertices  or  one  or  more  axes  of  symmetry. 
General  Schwarz-ChristofTcl  maps  do  not  appear  to  have  been  used  as  a 
computational  tool,  although  experiments  have  been  made  in  computing 
them. 


3.  Numerical  computation  of  the  Schwarz-Ghristoffel  Transformation 

In  the  early  days  of  computers,  when  a number  of  relatively  pure 
mathematicians  were  growing  interested  in  computational  mathematics, 
the  numerical  computation  of  conformal  maps  in  general  and  Schwarz- 
ChristolTcl  transformations  in  particular  received  a flurry  of  attention.  As 
early  as  1949,  the  National  Bureau  of  Standards  sponsored  a symposium  on 
numerical  conformal  mapping.  It  was  too  early,  however,  for  algorithms  to 
result  from  this  period  which  we  could  now  consider  practical. 

In  more  recent  years  interest  in  numerical  conformal  mapping  has  been 
modest.  Gaier  [1964]  produced  a comprehensive  work  describing  methods 
for  various  problems  in  constructive  conformal  mapping.  For  the  Schwarz- 
ChristofTcl  problem,  he  proposed  determining  the  accessory  parameters  Zk 
by  setting  up  a constrained  nonlinear  system  of  N — 3 equations  relating 
(1.4)  to  the  known  distances  |ufc — tuy|,  and  solving  it  iteratively  by  Newton’6 
method  [Gaier,  p.  171].  Such  a procedure  has  been  tried  by  at  least  three  6ets 
of  people:  [Meyer,  1979],  [Howe,  1973],  and  [Vecheslavov&Kokoulin,1973]. 

The  present  work  follows  Gaier  and  others  in  formulating  the  parameter 
problem  as  a constrained  nonlinear  system  of  equations.  We  believe  that 
this  is  the  first  fully  practical  program  for  computing  Schwarz-ChristofTel 
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transformations,  however,  and  the  first  which  is  capable  of  high  accuracy 
without  exorbitant  cost. 

One  innovation  which  makes  accurate  but  cheap  computations  possible 
here  is  the  use  of  a compound  form  of  Gauss-Jacobi  quadrature  to  evaluate 
the  integral  in  (1.4).  The  evaluation  of  this  integral  is  central  in  all  Schwarz- 
ChristolTcl  computations,  both  in  determining  the  accessory  parameters  and 
in  evaluating  the  map  and  its  inverse  once  the  accessory  parameters  are 
known.  We  have  found  that  a straightforward  application  of  Gau66-Jacobi 
quadrature,  as  6ome  others  have  U6ed,  can  achieve  only  very  low  accuracy 
in  realistic  problems,  and  we  have  developed  a compound  form  of  Gauss- 
Jacobi  quadrature  to  get  around  this  difficulty  (see  11.3). 

A second  innovation  here  is  that  the  computation  may  be  performed 
not  just  for  bounded  polygons,  but  for  polygons  with  any  number  of  vertices 
at  infinity.  This  is  made  possible  by  taking  the  unit  disk  as  the  model 
domain  rather  than  the  upper  half  plane,  which  others  have  used,  and 
evaluating  complex  contour  integrals  within  the  disk  rather  than  only  along 
the  boundary.  The  ability  to  handle  unbounded  polygons  is  important  for 
applications,  since  one  of  the  attractions  of  conformal  mapping  is  that  it 
can  reduce  an  unbounded  problem  domain  to  a bounded  one. 

The  treatment  of  the  constraints  in  the  nonlinear  system  is  a third 
new  feature  in  this  work.  We  have  employed  a simple  change  of  variables 
to  eliminate  these  constraints  directly.  This  approach  appears  to  be  more 
efficient  than  other  techniques  which  have  been  tried  (see  [Howe,  1973]  and 
[Vcchcslavov&Kokoulin,1973]),  and  eliminates  the  need  for  an  initial  guess 
of  the  accessory  parameters. 

We  have  depended  in  several  places  on  the  U6e  of  a sophisticated  library 
of  “black  box"  numerical  routines.  Library  programs  come  into  play  here 
for  Gauss-Jacobi  quadrature,  for  the  solution  of  the  nonlinear  system,  and 
for  the  solution  of  an  ordinary  differential  equation.  Others  have  been  used 
in  various  experiments  with  applications.  The  Schwarz-Christoffel  problem 
is  essentially  a simple  problem  numerically  once  the  machinery  i6  in  place, 
but  it  is  only  in  recent  years  that  this  kind  of  numerical  machinery  has 
begun  to  be  broadly  available. 
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II.  DETERMINATION  OF  THE  ACCESSORY  PARAMETERS 


1.  Formulation  a « a constrained  nonlinear  system  (subroutine  SCFUN) 

The  first  matter  to  be  settled  in  formulating  the  parameter  problem 
numerically  is,  what  parameters  in  the  map  (1.4)  6hall  we  fix  at  the  outset 
to  determine  the  Schwarz-Christoffel  transformation  uniquely?  One  choice 
would  be  to  fix  three  of  the  boundary  points  s^:  say,  *i  “ 1,  *2  ■=  *,*N  “ — »• 
This  normalization  has  the  advantage  that  the  resulting  nonlinear  system 
has  size  only  ( N — 3)-by-(N  — 3),  which  for  a typical  problem  with  N*=8 
may  lead  to  a solution  in  less  than  half  the  time  that  a method  involving 
an  (N — l)-by-(N  — I)  system  requires.  Nevertheless,  we  have  chosen  here 
to  normalize  by  the  conditions: 

***=!  (2.1) 

we  = arbitrary  point  within  P 

which  lead  to  an  ( N — 1)— by— (Af  — 1)  system.  This  choice  i6  motivated 
by  considerations  of  numerical  scaling:  it  allows  the  vertices  to  distribute 
themselves  more  evenly  around  the  unit  circle  than  they  might  otherwise. 
(An  earlier  version  of  the  program  mapped  from  the  upper  half  plane  instead 
of  the  unit  disk,  but  was  rejected:  once  points  **  began  appearing  far  from 
the  origin  at  x *=  10\  scaling  became  a problem.)  After  a map  has  been 
computed  according  to  any  normalization,  it  is  of  course  an  easy  matter  to 
transform  it  analytically  to  a different  domain  or  a different  normalization 
by  a Mobius  transformation. 

Now  the  nonlinear  system  must  be  formulated.  The  final  map  must 
satisfy  N complex  conditions, 

, 1 < k < N.  (2.2) 

These  amount  to  2 N real  conditions  to  be  satisfied,  but  they  are  heavily  over- 
determined,  for  the  form  of  the  Schwarz-Christoffel  formula  (1.4)  guarantees 
that  the  angles  will  be  correct  no  matter  what  accessory  parameters  are 
chosen.  We  must  reduce  the  number  of  operative  equations  to  N — 1.  This 
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i6  a tricky  matter  when  unbounded  polygons  are  allowed,  for  one  must  be 
careful  that  enough  information  about  the  polygon  P is  retained  that  no 
degrees  of  freedom  remain  in  the  computed  solution. 

We  proceed  a6  follows.  First,  we  require  that  every  connected  com- 
ponent of  P contain  at  least  one  vertex  to*.  Thus  even  an  infinite  straight 
boundary  must  be  considered  to  contain  a (degenerate)  vertex.  This  restric- 
tion eliminates  any  translational  degrees  of  freedom.  Second,  at  least  one 
component  of  P must  in  fact  contain  two  finite  vertices,  and  wn  and  tuj  will 
be  taken  to  be  two  such.  Th»  restriction  eliminates  rotational  degrees  of 
freedom. 

Now  define 


n(*  - ~) 


(2.3) 


where  xn  = 1 is  fixed  permanently  by  (2.1).  Next,  impose  the  complex 
condition  (real  equations  1,2) 


iuj  — wc 


(2.4 a) 


This  amounts  to  two  real  equations  to  be  satisfied. 

Denote  by  r*i , ...,  Tm  the  distinct  connected  components  of  P,  numbered 
in  counterclockwise  order.  For  each  t > 2,  impose  one  more  complex  con- 
dition: if  Zkf  is  the  last  vertex  of  in  the  counterclockwise  direction,  then 
(real  equations  3,4,. ..,2m) 


Wk,  — wc  = 


(2  46) 


Finally,  N — 2m  — 1 conditions  of  side  length  are  imposed.  For  each 
pair  (zk,Zk+i)  beginning  at  k = 1 and  moving  counterclockwise,  where  both 
vertices  are  finite,  we  require  (real  equations  2m -f- 1 — 1) 

. /•**+!  N / J 

i“»+* - «*i  - p i n(* - ^ 


\-fii  I 


(2-4c) 
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until  a total  of  N — 1 conditions  have  been  imposed.  If  P contains  at  least 
one  vertex  at  infinity,  then  every  bounded  side  will  have  been  represented  in 
a condition  of  the  form  (2.4c)  except  for  the  side  (wn,  u>i),  which  is  already 
taken  care  of  by  (2.1)  and  (2.4a).  If  P is  bounded,  then  the  last  two  sides 
in  counterclockwise  order — (u>A/_a»  to^_i)  and  (tu^— i,  w^) — will  not  be  so 
represented. 

We  have  not  stated  over  what  contours  the  integrals  of  eqs.  (2.4)  are 
defined.  This  does  not  matter  mathematically,  as  the  integrand  is  analytic, 
but  it  may  matter  numerically.  In  this  work  we  have  evaluated  them  always 
over  the  straight  line  segment  between  the  two  endpoints,  a procedure  which 
poses  no  domain  problems  6ince  the  unit  disk  is  strictly  convex.  Figure  2.1 
illustrates  what  contours  are  involved  in  computing  the  integrals  in  (2.3) 
and  (2.4),  for  a sample  case  with  N = 10,  m = 3. 

The  nonlinear  system  is  now  determined,  and  its  unique  solution  will 
give  the  unknown  parameters  C and  z\,...,zn—\  for  theSchwars-Christoffel 
mapping.  We  must,  however,  take  notice  of  two  special  cases  in  which  the 
solution  is  not  completely  determined  by  eqs.  (2.4).  It  was  remarked  that 
if  P is  bounded,  then  nowhere  in  eqs.  (2.4)  does  the  point  tu/v—  i appear.  If 
/?N_i  jk  — 1 or  0,  then  this  omission  is  of  no  consequence,  for  the  geometry 
of  the  problem  forces  to  be  correct.  If  Pn—\  = 0 or  — 1,  however,  then 
iu/v— i is  not  determined  a priori.  The  former  case  is  of  little  consequence, 
for  since  =ss  0 the  value  taken  fora^v—i  ha6  no  effect  on  the  computed 
mapping,  as  may  be  seen  in  (1.4),  nor  is  there  any  purpose  in  including  to/v— i 
among  the  vertices  of  P in  the  first  place.  (Still,  there  may  be  problems 
in  solving  the  system  (2.4)  numerically,  tor  it  i6  now  underdetermined.) 
The  latter  case,  Ps—i  = — 1,  is  more  serious,  and  must  be  avoided  in  the 
numbering  of  the  vertices  it*. 


2.  Transformation  to  an  unconstrained  system  (subroutine  YZTRAN) 

The  nonlinear  system  (2.4)  ostensibly  involves  N — 1 complex  unknown 

points  si i on  the  unit  circle.  In  dealing  with  such  a system,  we 

naturally  begin  by  considering  not  the  points  themselves,  but  their  argu- 
ments Ok,  given  by 

0<«»<2*.  (2.5) 

Now  the  system  depends  on  N — 1 real  unknowns,  and  the  solution  in  terms 
of  the  Ok  is  fully  determined. 

However,  the  system  (2.4)  as  it  stands  must  be  subject  to  a set  of  strict 
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Figure  2.1  - Contours  of  integration  within  the  disk.  A sample  Schwars- 
Christoffcl  problem  is  shown  with  N 10  vertices  of  which  m = 3 vertices  are 
at  infinity,  illustrating  what  integrals  are  computed  to  evaluate  the  system  (2.4): 


1 radial  integral  along  (0  — Z\q)  defines  C (eq.  2.3) 


1 radial  integral  along  (0  — x{]  determines  two  real  equations  to  fix  u>i  (eq. 
2.4a) 


• 2 radial  integrals  along  (0  — 25)  and  (0  — 27)  determine  four  real  equations 
to  fix  tt*5  and  to?  (eq.  2.4b) 


• 3 chordal  integrals  along  (zj  — Z4),  (24  — zg),  and  (*g — *jq)  determine  three 
real  equations  to  fix  |tM<  — u^lJu*  — W4J,  and  |u^o  — u*|  (eq.  2.4c) 


TOTAL:  N — 1 = 9 real  equations 
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inequality  constraints, 

0 < fl/b  < 0fc+i » (2.6) 

which  embody  the  fact  that  the  vertices  ifc  must  lie  in  ascending  order  coun- 
terclockwise around  the  unit  circle.  To  solve  the  system  numerically,  it  is 
desirable  to  eliminate  these  constraints  somehow.  We  do  this  by  transform- 
ing eqs.  (2.4)  to  a system  in  N — 1 variables  yi,...,yjv— i,  defined  by  the 
formula 


Vi-'oih^Z'  (s-3 * * * 7) 

where  Oq  and  On,  two  different  names  for  the  argument  of  mn  1,  are  taken 

for  convenience  as  0 and  2*,  respectively. 

At  each  iterative  step  in  the  solution  of  the  nonlinear  system  (2.4), 
we  begin  by  computing  a set  of  angles  {0*}  and  then  vertices  {z*}  from 
the  current  trial  set  {]/*}.  This  is  easy  to  do,  though  not  immediate  since 
the  equations  (2.7)  are  coupled.  In  this  way  the  problem  is  reduced  to  one 
of  solving  an  unconstrained  nonlinear  system  of  equations  in  N — 1 real 
variables. 


3.  Integration  by  compound  Gauss-Jacobi  quadrature  (subroutine  ZQUAD) 

The  central  computation  in  solving  the  parameter  problem,  and  indeed 

in  all  Schwarz-ChristofTel  computations,  is  the  numerical  evaluation  of  the 
Schwarz-Christoffel  integral  (1.4)  along  some  path  of  integration.  Typically 
one  or  both  endpoints  of  this  path  are  prevertices  ifc  on  the  unit  circle,  and 
in  this  case  a singularity  of  the  form  (1  — is  present  in  the  integrand 
at  one  or  both  endpoints. 

A natural  way  to  compute  such  integrals  quickly  is  by  means  of  Gauss- 
Jacobi  quadrature  (see  pavis  Sc  Rabinowitz,  1975],  p.  75).  A Gauss-Jacobi 
quadrature  formula  is  a sum  «>»/(*•),  where  the  weights  u*  and 

nodes  x,  have  been  chosen  in  such  a way  that  the  formula  computes  the 
integral  /(x)(  1 — x)°(l  + x)fldx  exactly  for  /(x)  a polynomial  of  as 
high  a degree  as  possible.  Thus  Gauss-Jacobi  quadrature  is  a generalization 
of  pure  Gaussian  quadrature  to  the  case  where  singularities  of  the  general 
form  (1  — x)°(l  + *)*  ( a,{J  > —1)  are  present.  The  required  nodes  and 
weights  can  be  computed  numerically;  we  have  U6ed  the  program  GAUSSQ 
by  Golub  and  Welsch  [Golub  Sc  Welsch,1969]  for  this  purpose. 
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Gauss-  Jacobi  quadrature  appears  made-to-order  for  the  Schwarz-Chris- 
tofTel  problem,  and  at  least  two  previous  experimenters  have  used  it  or  or 
a closely  related  technique  ([Howe, 1973],  [Vecheslavov  & Kokoulin,1973]). 
We  began  by  doing  the  6ame,  and  got  good  results  for  many  polygons  with 
a small  number  of  vertices.  In  general,  however,  we  found  this  method  of 
integration  very  inaccurate.  For  a typical  sample  problem  with  N «=  12 
and  NPTS  — 10,  it  produced  integrals  accurate  to  only  about  10  2,  and 
it  docs  much  worse  if  one  chooses  polygons  designed  to  be  troublesome. 

What  goes  wrong  is  a matter  of  resolution.  Consider  a problem  like  the 
one  show'n  in  Figure  2.2.  We  wish  to  compute  the  integral  (1.4)  along  the 
Bcgmcnt  from  z±  to  some  point  p.  (In  the  parameter  problem  p might  be  0 or 
Zk—\\  in  later  computations  it  might  be  any  point  in  the  disk.)  Now  direct 
application  of  a Gauss-Jacobi  formula  will  involve  sampling  the  integrand 
at  only  NPTS  nodes  between  zk  and  p.  If  the  singularity  z^i  is  so  close  to 
the  path  of  integration  that  the  distance  t = |*fc+i  — is  comparable  to 
the  distance  between  nodes,  then  obviously  the  Gauss-Jacobi  formula  will 
yield  a very  poor  result.  It  turns  out  that  in  Schwarz-Chri6toffel  problems 
the  correct  spacing  of  prevertices  z*  around  the  unit  circle  i6  typically  very 
irregular,  so  the  appearance  of  this  problem  of  resolution  is  the  rule,  not 
the  exception.  (See  examples  in  V.) 

To  maintain  high  accuracy  without  giving  up  much  speed,  we  have 
switched  to  a kind  of  compound  Gauss-Jacobi  quadrature  (see  [Davis  Sc 
Rabinowitz,  1975],  p.  56).  We  adopt,  somewhat  arbitrarily,  the  following 
quadrature  principle: 


No  singularity  zk  shall  lie  closer  to  an  interval  of 
integration  than  half  the  length  of  that  interval. 


To  achieve  this  goal,  the  quadrature  subroutine  ZQUAD  must  be  able  to 
divide  an  interval  of  integration  into  shorter  subintervals  as  necessary,  work- 
ing from  the  endpoints  in.  On  the  short  subinterval  adjacent  to  the  endpoint 
Gauss-Jacobi  quadrature  will  be  applied;  on  the  longer  interval  (or  intervals) 
away  from  the  endpoint  pure  Gaussian  quadrature  will  be  applied.  The 
eiTect  of  this  procedure  is  that  number  of  integrand  evaluations  required  to 
achieve  a given  accuracy  i6  reduced  from  0(£)  to  0(log2  £). 

Figure  2.2  shows  the  intervals  of  integration  that  come  into  play  in 
compound  Gauss-Jacobi  quadrature.  For  a plot  comparing  the  accuracy  of 
simple  and  compound  Gauss-Jacobi  quadrature  in  another  typical  problem, 
see  IV. 1. 
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Figure  2.2  - Compound  Gauss-Jacobi  quadrature.  Division  of  an  interval 
of  integration  into  eubintcrvals  to  maintain  desired  resolution. 


With  the  use  of  compound  Gauss-Jacobi  quadrature,  we  now  achieve 
high  accuracy  in  little  more  than  the  time  that  direct  Gauss-Jacobi  quad- 
rature takes.  This  i6  possible  because  only  a minority  of  integrals  have  a 
singularity  close  enough  that  subdivision  of  the  interval  of  integration  is  re- 
quired. In  the  12-vertcx  example  mentioned  above,  the  switch  to  compound 
Gauss-Jacobi  integration  decreased  the  error  from  10~ 3 to  2 • 10— 7. 

There  remains  one  circumstance  in  which  integration  by  compound 
Gauss-Jacobi  quadrature  as  described  here  is  unsuccessful.  This  is  the  case 
of  an  integration  interval  with  one  endpoint  quite  near  to  some  prevertex 
Zk  corresponding  to  a vertex  w*  = oo.  We  cannot  evaluate  such  an  integral 
by  considering  an  interval  which  begins  at  Zk,  for  the  integral  would  then 
be  infinite.  The  proper  approach  to  this  problem  is  probably  the  use  of 
integration  by  part6,  which  can  reduce  the  singular  integrand  to  one  that 
is  not  infinite.  Depending  on  the  angle  one  to  three  applications  of  in- 
tegration by  parts  will  be  needed  to  achieve  this.  We  have  not  implemented 
this  procedure. 

The  subtlety  of  the  integration  problem  in  Schwarz-ChristofTel  com- 
putations i6  worth  emphasizing.  It  is  customary  to  dispatch  the  integration 
problem  as  quickly  as  possible,  in  order  to  concentrate  on  the  "difficult"  ques- 
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t ion b:  computation  of  accessory  parameters  and  inversion  of  the  Schwarz- 
ChristofTcl  map.  We  believe,  however,  that  the  more  primary  problem  of 
computing  Schwarz-Christoffel  integrals — the  "forward"  problem — should 
always  remain  a central  concern.  Any  numerical  approach  to  the  parameter 
problem  or  the  inversion  problem  is  likely  to  employ  an  iterative  scheme 
which  depends  at  each  6tep  on  an  evaluation  of  the  integral  (1.4),  and  so 
the  results  can  only  be  as  accurate  as  that  evaluation. 


4.  Solution  of  system  by  packaged  solver  (subroutine  SGSOLV) 

The  unconstrained  nonlinear  system  is  now  in  place  and  ready  to  be 
solved.  For  this  purpose  we  employ  a library  subroutine:  NS01A,  by  M.J.D. 
Powell  ([Powell,  1968]),  which  uses  a steepest  descent  search  in  early  itera- 
tions if  necessary  followed  by  a variant  of  Newton's  method  later  on.  (The 
routine  does  not  use  analytic  derivatives.)  It  is  assumed  that  a variety  of 
other  routines  would  have  served  comparably  well. 

We  make  no  attempt  to  tailor  the  numerical  solution  procedure  to  the 
particular  Schwarz-Christoffel  problem  under  consideration.  In  particular, 
nil  iterations  begin  with  the  trivial  initial  estimate  y*  = 0 (1  < fc  < N — 

1).  This  corresponds  to  trial  vertices  spaced  evenly  around  the  unit  circle. 
The  following  input  parameters  to  NS01A  have  generally  remained  fixed: 
DSTEPe=10~ 8 (step  size  used  to  estimate  derivatives  by  finite  differences), 
DMAX  = 10  (maximum  step  size),  MAXFUN  «=  15 [N  — 1)  (maximum 
number  of  iterations). 

A fourth  parameter,  EPS,  defines  the  convergence  criterion — how  large 
a function  vector  (square  root  of  6um  of  squares  of  functions  values)  will 
be  considered  to  be  satisfactorily  close  to  zero.  We  have  most  often  taken 
10— 8 or  10~ 14  here.  The  choice  of  EPS  is  not  very  critical,  however,  as 
convergence  in  NS01A  is  generally  quite  fast  in  the  later  stages. 

In  the  course  of  thi6  work  about  a hundred  Schwarz-Christoffel  trans- 
formations have  been  computed,  ranging  in  complexity  from  N = 3 to 
N = 18.  NS01A  has  converged  successfully  to  an  accurate  solution  in  all 
of  these  trials.  Section  V.l  gives  a series  of  plots  showing  this  convergence 
graphically  for  a simple  example. 
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III.  COMPUTATION  OF  THE  S-C  MAP  AND  ITS  INVERSE 


Determining  the  accessory  parameters  is  the  most  formidable  task 
in  computing  numerical  Schwarz-Christoffel  transformations.  Once  this  U 
done,  evaluation  of  the  map  and  of  its  inverse  follow  relatively  easily.  The 
foundation  of  these  computations  continues  to  be  compound  Gauss>Jacobi 
quadrature. 


1.  From  disk  to  polygon:  w = w[z)  (subroutine  WSC) 

To  evaluate  the  forward  map  ti>(s)  for  a given  point  x in  the  disk  or  on 
the  circle,  we  must  compute  the  integral 

u.-^+cf  n(i-- 

Ja°  y=i ' 

with  uaq  = iu(ao),  where  the  endpoint  zo  may  be  any  point  in  the  closed  disk 
at  which  the  image  w[zo)  is  known  and  not  infinite.  Three  possible  choices 
for  zo  suggest  themselves — 

(1)  zo  = 0;  hence  wo  = we‘, 

(2)  zo  = Zk  for  some  k ; hence  u*  = ufr,  a vertex  of  P; 

(3)  zo  — some  other  point  in  the  disk  at  which  w has  previously  been 
computed. 

In  cases  (1)  and  (3),  neither  endpoint  has  a singularity,  and  an  evaluation  of 
(3.1)  by  compound  Gauss-Jacobi  quadrature  reduces  to  the  U6e  of  compound 
Gauss  quadrature.  In  case  (2)  a singularity  of  the  form  (1  — z/zk)~p  is 
present  at  one  of  the  endpoints  and  the  other  endpoint  has  no  singularity. 
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The  best  rule  for  computing  tt>(z)  is:  if  x is  close  to  a singular  point  a* 
(but  not  one  with  = oo),  use  method  (2);  otherwise,  use  method  (1).  In 
either  case  we  employ  compound  Gauss-Jacobi  quadrature,  taking  normally 
the  same  number  of  nodes  as  was  used  in  solving  the  parameter  problem. 
By  this  procedure  we  evaluate  w(z)  readily  to  "full"  accuracy — that  is,  the 
accuracy  to  which  the  accessory  parameters  have  been  computed,  which  is 
directly  related  to  the  number  of  points  chosen  for  Gauss-Jacobi  quadrature 
(sec  IV.  1).  Quadrature  nodes  and  weights  need  only  be  computed  once,  of 
course. 

We  should  emphasize  that  even  in  the  vicinity  of  a singularity  z/t,  the 
evaluation  of  the  map  u>  = u>(z)  is  inherently  very  accurate.  Thi6  very 
satisfactory  treatment  of  singular  vertices  is  a considerable  attraction  of 
the  Schwarz-Christoffel  approach  for  solving . problems  of  Laplace  type. 
In  particular,  in  a potential  problem  the  Schwarz-Christoffel  transforma- 
tion "automatically"  handles  the  singularities  correctly  at  any  number  of 
reentrant  corners. 


2.  From  polygon  to  disk:  z = z(w)  (subroutine  ZSC) 

For  computing  the  inverse  mapping  z = z{w)  at  least  two  possibilities 
exist,  both  of  them  quite  powerful.  The  most  straightforward  approach  is 
to  view  the  formula  w[z)  = w as  a nonlinear  equation  to  be  solved  for  z, 
given  6omc  fixed  value  w.  The  solution  may  then  be  found  iteratively  by 
Newton's  method  or  a related  device.  w[z)  should  be  evaluated  at  each  step 
of  such  a process  by  compound  Gauss-Jacobi  quadrature  along  a straight 
line  segment  whose  initial  point  remains  fixed  throughout  the  iteration. 

An  alternative  approach  is  to  invert  the  Schwarz-Christoffel  formula, 


to  yield  the  formula 


fi-is-sr 


This  inversion  is  possible  because  iu  «=  w(z)  is  a conformal  mapping,  which 
means  \dw/dz\  > 0 everywhere.  (3.2)  may  now  be  thought  of  a6  an  ordinary 


differential  equation  (o.d.e.), 


(3.3) 

in  one  complex  variable  w.  If  a pair  of  values  (ao,  wo)  is  known  and  the  new 
value  z = z[w)  is  sought,  then  z may  be  computed  by  applying  a numerical 
o.d.e.  solver  to  the  problem  (3.3),  taking  as  a path  of  integration  any  curve 
from  itb  to  w which  lies  within  the  polygon  P. 

In  our  program  we  have  chosen  to  combine  these  two  methods,  using  the 
second  method  to  generate  an  initial  estimate  for  use  in  the  first.  We  begin 
with  the  o.d.e.  formulation,  using  the  code  ODE  by  Shampine  and  Gordon, 
and  for  convenience  we  integrate  whenever  possible  along  the  straight  line 
segment  from  wc  to  w.  (ODE,  like  most  o.d.e.  codes,  is  written  for  problems 
in  real  arithmetic,  so  that  we  mu6t  first  express  (3.2)  as  a system  of  first- 
order  o.d.e.'s  in  two  real  variables.)  Since  P may  not  be  convex,  more  than 
one  line  segment  step  may  be  required  to  get  from  wo  to  w in  this  way.  It 
will  not  do  to  take  tuo  = w*  for  some  vertex  u*  without  special  care,  because 
(3.2)  is  singular  at  u>*. 

From  ODE  we  get  a rough  estimate  z of  z(w),  accurate  to  roughly  10— 2. 
This  estimate  is  now  used  as  an  initial  guess  in  a Newton  iteration  to  solve 
the  equation  w(z)  = w.  This  method  is  faster  than  the  o.d.e.  formulation  for 
getting  a high-accuracy  answer.  More  important,  it  is  based  on  the  central 
Gauss-Jacobi  quadrature  routine,  unlike  the  o.d.e.  computation. 

In  summary,  wc  compute  the  inverse  map  z = z(w)  rapidly  to  full 
accuracy  by  the  following  steps: 

(1)  Solve  (3.2)  to  low  accuracy  with  package  ODE,  integrating  when- 
ever possible  along  the  line  segment  from  wc  to  w ; call  the  result 

(2)  Solve  the  equation  w{z)  « u>  for  z by  Newton's  method,  using  z 
as  an  initial  guess. 
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IV.  ACCURACY  AND  SPEED 


1.  Accuracy 

The  central  computational  step  is  the  evaluation  of  the  Schwarz-Chris- 
tofTel  integral,  and  the  accuracy  of  this  evaluation  normally  determines 
the  accuracy  of  the  overall  computation.  As  a consequence  of  the  quadra- 
ture principle  adopted  in  II. 3 — that  no  quadrature  interval  shall  be  longer 
than  twice  the  distance  to  the  nearest  singularity  z* — the  compound  Gauss- 
Jacobi  formulation  achieves  essentially  the  full  accuracy  typical  of  Gaussian 
quadrature  rules  operating  upon  smooth  integrands.  That  is,  the  number  of 
digits  of  accuracy  i6  closely  proportional  to  NPTS,  the  number  of  quadrature 
nodes  per  half-interval,  with  a very  satisfactory  porportionality  constant  in 
practice  of  approximately  I. 

It  is  important  not  only  to  be  capable  of  high  accuracy,  but  to  be 
able  to  measure  how  much  accuracy  one  has  in  fact  achieved  in  a given 
computation.  To  do  this  we  employ  a subroutine  TEST,  which  is  regularly 
called  immediately  after  the  parameter  problem  is  solved.  Given  a computed 
set  of  accessory  parameters  C and  {**}  , TEST  computes  the  distances 
|u*  — it»c|  for  each  u;*  ^ oo  and  the  distances  \wk—i  — u>*+i|  for  each 
ufc  = oo,  making  use  of  the  standard  subroutine  ZQUAD  for  compound 
Gauss-Jacobi  quadrature.  The  numbers  obtained  are  compared  with  the 
exact  distances  specified  by  the  geometry  of  the  polygon,  and  the  maximum 
error,  RADEMX,  is  printed  as  an  indication  of  the  magnitude  of  errors  in 
the  converged  solution.  It  is  now  probable  that  subsequent  computations  of 
w(z)  or  z( w)  will  have  errors  no  greater  than  roughly  RADEMX. 

Most  often  wc  have  chosen  to  use  an  8-point  quadrature  formula.  Since 
each  interval  of  integration  is  initially  divided  in  half  by  subroutine  ZQUAD, 
this  means  in  reality  at  least  16  nodes  per  integration.  With  this  choice 
RADEMX  consistently  has  magnitude  ~10— ; 8 for  polygons  on  the  scale  of 
unity. 

Figure  4.1  gives  an  indication  of  the  relationship  between  number  of 
quadrature  nodes  and  error  RADEMX;  it  shows  RADEMX  as  a function  of 
NPTS  for  a 6-gon  which  is  shown  at  the  top  of  the  next  page.  Two  curves 
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are  shown:  one  for  simple  Gau66-Jacobi  quadrature,  and  one  for  compound 
Gauss-Jacobi  quadrature.  The  exact  quantities  here  should  not  be  taken  too 
seriously;  examples  could  easily  have  been  devised  to  make  the  difference  in 
performance  of  the  two  quadrature  methods  much  smaller  or  much  greater. 


2.  Speed 

Any  application  of  Schwarz-Christoflfel  transformations  consists  of  a 
sequence  of  steps: 

INIT  - set  up  problem 

Q1NIT  - compute  quadrature  nodes  and  weights 

SCSOLV  - solve  parameter  problem 

TEST  - estimate  accuracy  of  solution 

ZSC,  WSC,  etc.  - compute  forward  and  inverse  transformations  in 
various  applications 

Among  these  tasks  INIT,  QINIT,  and  TEST  all  take  negligible  amounts 
of  time  relative  to  the  other  computations:  typically  less  than  0.1  secs,  on 
the  IBM  370/168  for  INIT  and  QINIT,  and  for  TEST  a variable  time  that 
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is  usually  less  than  5%  of  the  time  required  by  SCSOLV.  What  remains  are 
three  main  time  consumers:  SCSOLV,  ZSC,  and  WSC. 

We  begin  with  WSC,  which  performs  the  central  evaluation  of  (1.4) 
by  compound  Gauss- Jacobi  quadrature.  This  evaluation  takes  time  propor- 
tional to  NPTS  (the  number  of  quadrature  nodes)  and  to  N (the  number  of 
vertices).  The  first  proportionality  i6  obvious,  and  the  second  results  from 
the  fact  that  the  integrand  of  (1.4)  U an  JV-fold  product.  Very  roughly,  we 
may  estimate 

time  to  solve  w = w(z) : 0.25  • NPTS  • N msec.  (4.1a) 

for  double  precision  computations  on  the  IBM  370/168.  Taking  a typical 
value  of  NPTS=8,  which  normally  leads  to  8-digit  accuracy,  (4.1a)  may  be 
rewritten 

(4.16) 

For  the  minority  of  cases  in  which  the  interval  must  be  subdivided  to 
maintain  the  required  resolution,  these  figures  will  be  larger. 

To  estimate  the  time  required  to  solve  the  parameter  problem,  we  com- 
bine (4.1)  with  an  estimate  of  how  many  integrals  must  be  computed  in  the 
course  of  solving  this  problem.  To  begin  with,  at  each  iteration  about  N 
integrals  arc  required  by  NS01A  (the  exact  number  depends  on  the  number 
of  vertices  at  infinity).  On  top  of  this,  it  is  a fair  estimate  to  say  that  4 N 
iterations  will  be  required  by  NS01A  to  achieve  a high-accuracy  solution. 
We  arc  therefore  led  to  the  estimate 

time  to  solve  parameter  problem:  NPTS  • N3  msec.  (4.2 a) 

or,  taking  again  NPTS=8, 

time  to  solve  parameter  problem:  8 N3  msec.  (4.26) 

These  estimates  correspond  fairly  well  with  observed  computation  times 
for  the  parameter  problem:  two  problems  with  N = 5 and  N «=  18  may 
be  expected  to  take  about  1 and  50  seconds,  respectively.  It  is  clear  that 
computing  a Schwarz-ChristofTel  transformation  becomes  quite  a sizeable 
problem  for  polygons  with  more  than  ten  vertices.  In  particular,  6uch  com- 
putations are  much  too  timc-con6uming  for  it  to  be  practical  to  approximate 
a curved  domain  by  a polygon  with  a large  number  of  vertices. 
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Finally,  we  must  consider  the  time  taken  by  subroutine  ZSC  to  invert 
the  Schwarz-ChristofTel  map.  This  too  is  proportional  to  NPTS,  and  quite 
problem  dependent.  We  estimate  very  roughly: 

time  to  solve  s = z(w) : NPTS  • N msec.  (4.3a) 

or,  with  NPTS— 8, 


(4.36) 


Note  that  inverting  the  Schwarz-ChristofTel  map  is  only  about  four 
times  as  time-consuming  as  computing  it  in  the  forward  direction. 

In  practice,  computational  applications  will  vary  considerably  in  the 
use  they  make  of  a Schwarz-ChristofTel  transformation  once  the  parameter 
problem  is  solved.  If  only  a few  dozen  applications  of  ZSC  or  WSC  are 
required,  then  the  computational  time  for  solving  the  parameter  problem 
will  dominate.  If  thousands  of  such  computations  are  needed,  on  the  other 
hand,  then  the  parameter  problem  may  become  relatively  insignificant.  The 
latter  situation  i6  most  likely  to  hold  when  plotting  is  being  done,  or  when 
a high-accuracy  solution  in  the  model  domain  is  to  be  computed  by  means 
of  finite  differences. 

In  summary,  high  accuracy  is  cheap  in  Schwarz-ChristofTel  transfor- 
mations; what  consumes  time  is  solving  problems  involving  a large  number 
of  vertices. 


time  to  solve  z = z{w) : 8 N msec. 


V.  COMPUTED  EXAMPLES  AND  APPLICATIONS 
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1.  Iterative  process  for  a single  example 

Figure  5.1  shows  graphically  the  process  of  convergence  from  the  ini- 
tial estimate  in  an  example  involving  a 4-gon.  Routine  NS01A  begins  by 
evaluating  the  function  vector  (2.4)  at  the  initial  guess,  then  at  each  of 
N — 1 input  vectors  determined  by  perturbing  the  initial  guess  by  the  small 
quantity  DSTEP  in  each  component.  As  a result,  the  first  N pictures  always 
look  almost  alike,  which  is  why  the  series  shown  begins  at  NEVAL=4  rather 
than  NEVAL=1.  Each  plot  shows  the  current  image  polygon  together  with 
the  images  of  concentric  circles  in  the  unit  disk  (which  appear  as  "contours") 
and  the  images  of  radii  leading  from  the  center  of  the  disk  to  the  current 
preverticcs  z 

These  pictures  have  a beautiful  bonus  feature  about  them:  they  may 
be  interpreted  as  showing  not  only  the  image  polygon  but  simultaneously 
the  domain  disk,  including  the  prevertices  Zk  along  the  unit  circle.  To  6ee 
this,  look  at  one  of  the  inner  "contour"  curves,  one  which  is  apparently 
circular,  and  the  radii  within  it.  Since  tu  = w[z)  is  a conformal  map  within 
the  interior  of  the  disk,  the  radii  visible  in  this  circle  must  intersect  at  the 
same  angles  as  their  preimages  in  the  domain  disk.  Thus  the  inner  part  of 
any  one  of  these  image  plots  is  a faithful  representation  on  a small  6cale  of 
the  circular  domain.  We  see  in  Figure  5.1  that  the  prevertices  are  equally 
spaced  around  the  unit  circle  initially  (NEVAL  «=  4),  but  move  rapidly  to 
a very  uneven  distribution.  This  behavior,  which  is  typical,  indicates  why 
the  use  of  a compound  form  of  Gauss-Jacobi  quadrature  is  60  important  (see 
n.3). 

The  sum-of-squares  error  in  solving  the  nonlinear  6y6tem  is  plotted  as 
a function  of  iteration  number  in  Figure  5.2,  for  the  same  4-vertex  example. 
Convergence  is  more  or  less  quadratic,  as  one  would  expect  for  Newton's 
method.  The  irregularity  at  iteration  19  is  caused  by  the  finite  difference 
step  size  of  10-8  used  to  estimate  derivatives,  and  would  have  been  repeated 
at  each  alternate  step  thereafter  if  the  iteration  had  not  terminated. 
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FUM-OF-SQUARES  ERROR 


— 


ITERATION  NUMBER 


Figure  5.2  - Rate  of  convergence.  Sum-of-iquarea  error  in  the  nonlinear 
ay  stem  (2.4)  aa  a function  of  iteration  number,  for  the  tame  problem 
aa  in  Figure  5.1. 
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2.  Sample  Schwars-Christoffel  maps 


Figures  5.3  and  5.4  show  plots  of  computed  Schwarz-Christoffel  maps 
for  representative  problems.  The  polygons  of  Figure  5.3  are  bounded  and 
those  of  Figure  5.4  are  unbounded.  Observe  that  contour  lines  bend  tightly 
around  reentrant  corners,  revealing  the  large  gradients  there,  while  avoiding 
the  backwater  regions  near  outward-directed  corners  and  vertices  at  infinity. 
Like  the  plots  of  Figure  5.1,  these  may  be  viewed  as  showing  simultaneously 
the  image  polygon  and  the  domain  disk. 

Figure  5.5  6how6  similar  plots  in  which  streamlines  rather  than  con- 
tour lines  have  been  plotted,  60  that  the  configuration  may  be  thought  of 
as  portraying  ideal  irrotational  fluid  flow  through  a two-dimensional  chan- 
nel. To  plot  these  streamlines  an  analytic  transformation  of  the  disk  to  an 
infinite  channel  with  straight  parallel  sides  was  used  in  conjunction  with  the 
Schwarz-ChristoiTel  transformation  from  the  disk  to  the  problem  domain. 


3.  Laplace’s  equation 

Conformal  maps  do  not  6olve  problems,  but  they  may  reduce  hard 
problems  to  easier  ones.  How  much  work  must  be  done  to  solve  the  easier 
problem  will  vary  considerably  with  the  application. 

(1)  In  the  best  of  circumstances,  the  original  problem  may  be  reduced 
to  a model  problem  whose  solution  is  known  exactly.  This  is  the 
case  in  the  fluid  flow  problems  of  Figure  5.5,  in  which  a crooked 
channel  may  be  mapped  to  an  infinite  straight  channel  of  constant 
width. 

(2)  If  a problem  of  Laplace’6  equation  with  pureDirichlct  or  Neumann 
boundary  conditions  can  be  mapped  conformally  to  a disk,  then 
Poisson’s  formula  or  Dini's  formula  [Kantorovich  & Krylov,  1958] 
provide  integral  representations  of  the  solution  at  each  interior 
point.  Such  integrals  may  be  evaluated  readily  on  the  computer 
to  yield  high  accuracy  solutions.  The  primary  disadvantage  of 
this  approach  is  that  a new  integral  mu6t  be  evaluated  for  each 
point  at  which  the  solution  is  desired. 

(3)  If  the  solution  will  be  required  at  many  points  in  the  domain, 
then  it  i6  probably  more  efficient  to  6olve  Laplace's  equation  by 
a trigonometric  expansion  of  the  form  ao  + rfc(a*  sin  £0  -j- 
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bit  cos  k0)\  coefficients  a*  and  6 * are  selected  so  as  to  fit  the  boundary 
conditions  closely.  A disadvantage  of  this  method  is  that  conver- 
gence of  the  expansion  may  be  slow  if  the  boundary  conditions 
are  not  smooth. 

(4)  Finally,  if  simpler  methods  fail,  a solution  in  the  model  domain 
may  be  found  by  a finite-difference  or  finite-element  technique. 
For  problems  of  Poisson's  equation  or  more  complicated  equations 
this  will  probably  normally  be  necessary. 

Figure  5.6  presents  an  example  of  type  (1).  We  are  given  an  infinite 
region  bounded  by  one  straight  boundary  fixed  at  potential  <p  — 0 and  one 
jagged  boundary  fixed  at  — 2.  We  may  think  of  this  as  an  electrostatics 
problem.  The  central  question  to  be  answered  computationally  will  be:  what 
are  the  voltage  and  the  electric  field  E ■=  — Vy?  at  a given  point,  either 
within  the  field  or  on  the  boundary? 

We  proceed  by  mapping  the  given  region  onto  the  disk  by  a Schwarz- 
Christoffcl  transformation,  then  analytically  onto  an  infinite  straight  chan- 
nel (as  in  the  examples  of  Figure  5.5).  In  the  straight  channel  and  E are 
known  trivially,  and  this  information  may  be  transferred  to  the  problem 
domain  through  a knowledge  of  the  conformal  map  that  connects  them  and 
of  its  (complex)  derivative.  We  omit  the  details,  which  are  straightforward. 

Figure  5.6b  shows  (FI  as  a function  of  x on  the  upper  and  lower  bound- 
aries of  the  region.  To  see  more  of  the  behavior  of  the  solution  field  near 
a reentrant  corner,  we  also  compute  the  field  at  three  points  near  3-j-  1.5s. 
These  results  are  given  in  Figure  5.6c. 


4.  Poisson’s  equation 

Consider  the  7-sided  region  shown  in  Figure  5.7a.  We  wish  to  solve 
Poisson's  equation 

A^(x,  y)  — ^ sin  2x(l  — 2(y  + l)2) 
on  this  region  subject  to  Dirichlet  conditions 

4>{*,  v)  — P{*,  V)  — ^ "in  2 x(y  + l)2 

on  the  boundary.  We  proceed  by  mapping  the  domain  to  the  disk  and 
solving  a transformed  problem  in  the  disk  in  polar  coordinates  by  means  of 
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(a)  Problem  domain:  region  between  two 
conducting  sheets 


(b)  Field  strength  along  the  top  boundary 

(solid  line)  and  bottom  boundary  (broken 
line) 
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(c)  Computed  potential  and  field  strength  at 
three  points  near  3 + 1.5i 


Figure  5.6 — Laplace  equation  example:  electric 
potential  and  field  between  two  infinite  sheets, 
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a second-order  fast  finite  difTerence  6olver  (PWSPLR,  by  P.  Swarztrauber 
and  R.  Sweet).  p[x,  y)  is  the  correct  solution  in  the  interior  as  well  as  on  the 
boundary,  so  we  can  determine  the  accuracy  of  the  numerical  solution. 

This  is  not  as  satisfactory  a procedure  as  was  available  for  Laplace 
equation  problems.  According  to  (1.2),  the  model  problem  here  is  Poisson’s 
equation  in  the  disk  with  an  altered  right  hand  side  containing  the  factor 
|/(z)|2,  where  / is  the  composite  map  from  the  disk  to  the  7-gon.  Two 
difficlutics  arise.  The  first  i6  that  to  set  up  the  transformed  equation  in  the 
disk,  p{wij)  must  be  computed  for  every  t = u >(z,y)  which  is  an  image  of 
a grid  point  in  the  disk.  This  is  time  consuming,  one  hundred  times  more 
so  in  this  experiment  than  the  fast  solution  of  Poisson’s  equation  once  it 
is  set  up.  Second,  |/f(z)| 2 is  singular  (unbounded,  in  this  example)  at  each 
prevertex  z^,  and  this  appears  to  interfere  with  the  second-order  accuracy 
which  we  would  like  to  observe.  The  table  in  Figure  5.7b  attests  to  both  of 
these  problems. 


5.  Eigenfrequencies  of  the  Laplace  operator 

Petter  Bjdrstad  (Computer  Science  Dept.,  Stanford  University)  has 
recently  combined  the  present  Schwarz-Christoffel  computation  with  a fast 
finite-difference  scheme  to  successfully  compute  eigenvalues  and  eigenvec- 
tors of  the  Laplacian  operator  on  polygonal  regions.  These  results  may  be 
interpreted  as  giving  the  normal  modes  and  frequencies  of  a thin  membrane 
in  two  dimensions,  or  of  a three-dimensional  waveguide  with  constant  cross- 
section.  This  work  will  be  reported  elsewhere. 


(a)  7-sided  problem  domain,  including  image  of  16x32 
finite-difference  arid  in  the  unit  disk 


Transformation 


Grid 

(rxfl) 

and  setup 
time 

Fast  Poisson 
solver  time 

Max.  error 

RMS 
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(b)  Computed  results  for  four  different  grids.  Time 
estimates  are  for  an  IBM  370/168. 


Figure  5.7—  Poisson  equation  example.  Problem  is 
transplanted  conformally  to  the  unit  disk  and  solved 
by  finite  differences. 
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A program  has  been  described  which  computes  accurate  Schwarz-Chris- 
toffel transformations  from  the  unit  disk  to  the  interior  of  a simply  connected 
polygon  in  the  complex  plane,  which  may  be  unbounded.  Key  features  of 
the  computation  have  been: 

(1)  Choice  of  the  unit  disk  rather  than  the  upper  half  plane  as  the 
model  domain,  for  better  numerical  scaling  (II.  1) 

(2)  Use  of  complex  contour  integrals  interior  to  the  model  domain 
rather  than  along  the  boundary,  making  possible  the  treatment 
of  unbounded  polygons  (II.  1) 

(3)  Use  of  compound  Gauss-Jacobi  quadrature  in  complex  arithmetic 
to  evaluate  the  Schwarz-Christoffel  integral  accurately  (n.3,HI.  1) 

(4)  Formulation  of  the  parameter  problem  as  a constrained  nonlinear 
system  in  N — 1 variables  (HI) 

(5)  Elimination  of  constraints  in  the  nonlinear  system  by  a simple 
variable  transformation  (13.2) 

(6)  Solution  of  the  system  by  a packaged  nonlinear  systems  solver; 
no  initial  estimate  required  (13.4) 

(7)  Computation  of  a reliable  estimate  of  the  accuracy  of  further 
computations,  once  the  parameter  problem  ha6  been  solved  (IV.  1) 

(8)  Accurate  evaluation  of  the  inverse  mapping  in  two  6teps  by  means 
of  a packaged  o.d.e.  solver  and  a packaged  complex  rootfinder 

(ni.2) 

Previous  efforts  at  computing  Schwarz-Christoffel  transformations  nu- 
merically include  [Cherednichenko  &,  Zhelankina,  1975],  [Hopkins  & Rob- 
erts, 1978],  [Howe,  1973],  [Meyer,  1979],  and  [Vecheslavov  & Kokoulin, 
1973).  The  present  work  differs  from  these  in  that  it  deals  directly  with 
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complex  arithmetic  throughout,  taking  the  unit  disk  rather  than  the  upper 
half  plane  as  the  model  domain  and  evaluating  complex  contour  integrals. 

This  makes  possible  the  computation  of  transformations  involving  general 
unbounded  polygons.  (Cherednichenko  & Zhelankina  [1975]  also  treat  un- 
bounded polygons,  by  a different  method.)  Two  other  important  differences 
are  the  use  of  compound  Gauss-Jacobi  quadrature,  and  the  application  of 
a change  of  variables  to  eliminate  constraints  in  the  nonlinear  system  ( (5), 
above).  We  believe  that  our  program  computes  Schwarz-Christoffel  trans- 
formations faster,  more  accurately,  and  for  a wider  range  of  problems  than 
previous  attempts. 

A variety  of  directions  for  further  work  suggest  themselves.  Here  are 
some  of  them: 

(1)  More  attention  should  be  paid  to  the  problem  of  inverting  the 
Schwarz-ChristofTel  map.  The  two-step  method  described  in  IQ. 2 
is  only  one  of  many  possibilities. 

(2)  The  program  could  easily  be  extended  to  construct  maps  onto  the 
exterior  of  a polygon — that  is,  the  interior  of  a polygon  whose 
interior  includes  the  point  at  infinity.  This  extension  would  be 
necessary  for  applications  to  airfoil  problems. 

(3)  It  should  not  be  too  great  a step  to  raise  the  present  program  to  the 
level  of  "software”  by  packaging  it  flexibly,  portably,  and  robustly 
enough  that  naive  users  could  apply  it  to  physical  problems. 

(4)  The  program  might  be  extended  to  handle  the  rounding  of  corners 
in  Schwarz-ChristofTel  transformations  [Henrici,1974].  What  about 
mapping  doubly  or  multiply  connected  polygonal  regions,  per- 
haps by  means  of  an  iterative  technique  which  computes  an  S-C 
transformation  at  each  step?  What  about  applying  S-C  transfor- 
mations to  eliminate  corners  in  the  conformal  mapping  of  curved 
domains? 


Most  important,  further  work  is  needed  in  the  direction  of  applications 
to  Laplace's  equation,  Poisson's  equation,  and  related  problems.  Irregular 
or  unbounded  domains  arc  generally  troublesome  to  deal  with  by  standard 
techniques,  particularly  when  singularities  in  the  form  of  reentrant  corners 
arc  present.  Schwarz-ChristofTel  transformations  oiler  a means  of  getting 
around  6uch  difficulties  in  a natural  way.  Much  more  experience  is  needed 
here. 


35 


APPENDIX:  PROGRAM  LISTING 


The  boundaries  of  this  program  are  not  sharply  defined,  for  the  configu- 
ration  changes  according  to  what  applications  are  being  treated.  The  present 
listing  includes  only  the  core  routines  used  to  6olve  the  parameter  problem 
and  to  evaluate  the  Schwarz-ChristofTel  function  and  its  inverse. 

An  experimental  copy  of  the  package  may  be  obtained  in  machine* 
readable  form  from  the  author. 

Control  program: 

SC 

Set-up: 

LNIT  initializes  variables  and  reads  input  data 
QINIT  computes  quadrature  nodes  and  weights 

Solution  of  parameter  problem: 

SCSOLV  controls  solution  of  parameter  problem 
YZTRAN  transforms  to  an  unconstrained  system 
SCFUN  sets  up  the  nonlinear  system  to  be  solved 
SCOUTP  prints  output  from  SCSOLV 
TEST  estimates  accuracy  of  computed  solution 

Compound  Gauss-Jacobi  quadrature: 

ZQUAD  divides  the  integral  into  two  halves 
ZQUAD1  evaluates  the  half-integral  (compound) 

DIST  finds  the  distance  to  the  nearest  singularity 
ZQSUM  sums  a Gauss-Jacobi  quadrature  rule 

Forward  and  inverse  S-C  map: 

WSC  evaluates  map  from  disk  to  polygon 
ZSC  evaluates  map  from  polygon  to  disk 
ZFODE  computes  initial  estimate 
ZNEWT  inverts  map  by  Newton's  method 

Miscellaneous  routines: 

ZPROD  evaluates  N-fold  Schwarz-ChristofTel  integrand 
FINITE  returns  "true"  if  the  argument  is  finite 
ENTER  begins  timing  of  the  current  subroutine 
EXIT  concludes  timing  of  the  current  subroutine 
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Library  routines  not  listed: 

GAUSSQ  (Golub&Welsch)  computes  Gauss-Jacobi  nodes  and  wts 
(called  by  QINIT) 

NS01A  (Powell)  solves  the  nonlinear  system 
(called  by  SCSOLV) 

ODE  (Shampine  Sc  Gordon)  solves  the  inverse  mapping  problem 
(called  by  ZSC) 


MAIN  PROGRAM 


PWCOPA!'  SC  - "SCHWARZ- CHRI STOP PEL" 


TH:j  PP'GPA.1  COMPUTES  THE  schwabz-christoppel  tiarspobhatior 
WHICH  S-RDS  THE  OMIT  DISH  TO  THE  INTERIOR  IP  THE  POLYGON 
Wm,...,W(H).  THIS  HAPPING  IS  OP  THE  PORN: 

Z N 

W * WC  ♦ C • IRT  PBOD  ( 1-Z/Z (K)  ) **BETAH  (R)  DZ  . 

0 K»1 


CD 


Z TO  CC  LYE  THE  PR  jBL  EH  WE  BEGIN  BY  PIRDIIG  THE  ACCESSOIY  PARAHETEBS 
Z — VERTTCPS  Z(K)  AND  CORSTART  C — POi  THE  HAP  OP  (1).  THIS  IS  DONE 
C BY  SOBprOTIRE  SCSOLP. 

C THE  IMAGE  POLYGON  HAY  BE  ORBOORDED;  PERMITTED  ANGLES  LIE  IN  THE 
C PANGE  -3. IE. 3ETAH (R) . LE. 1.  W (N)  ARD  W(1)  HOST  BE  FINITE. 

C WE  M-RHALIZ2  BY  THE  CONDITIONS: 

C Z(N)  * 1 (2.1) 

Z W(O.C)  « WC  (A  POINT  IN  THE  INTERIOR  OP  THE  POLYGON)  (2.2) 

C NOTATION: 

W(R)  - VERTEX  R OP  THE  IMAGE  POLYGON 
I 7 ( K)  - POINT  GR  THE  GRIT  DISK  HAPPED  TO  W (R) 

Z 3£T  AM  (R)  - NBGATIf E IP  EXTEPIOR  ANGLE  AT  i (R)  DIVIDED  BY  PI 

Z N - NHMPEB  OP  VERTICES  W (R) 

T MM  * F-1  - NUMBER  DP  UNKNOWN  POINTS:  Z ( 1)  , . . . , Z (N- 1) 

T MPTSQ  - NUMBER  Of  POINTS  POP  GAUSS-JACOBI  QUADRATURE 

Z Z'NP  - COMPLEX  I NPIRITY 

C LOCAL  OTIRES: 

C SC  - MAIN  PROGRAM 

C TNI T - INITIALIZES  CONSTANTS  AND  DEPIRES  PROBLEM 

C Or*!!  - COMPUTES  QUADRATURE  NODES  AND  WEIGHTS 
C SCSOLV  - COMPUTES  ACCESSORY  PARAMETERS  POR  S-C  HAP  (1) 

C YZ-PAN  - TPARSPORMS  UNKNOWNS  FROM  Y (R)  TO  Z(R) 

C 3CPUN  - NONLINEAR  SYSTEM  OP  EQUATIONS  TO  BE  SOLVED  BY  SCSOLV 

C SCOUTP  - PRINTS  OUTPUT  FFGH  SCSOLV 

C W ??  - COMPUTES  W(Z) 

C ZSC  - COMPUTES  Z (R) 

C PLTCOM  - DRAWS  PLOTS  CP  IMAGE  POLYGON  RITH  CONTOURS 
C ZPFOD  - COMPUTES  R-POLD  PRODUCT  IB  (1) 

Z 7QUA0  - SUMS  TO  EVALUATE  INTEGRAL  BY  GAUSS-JACOBI  QUADRATURE 

C PTIIIE  - RETURNS  TRUE  IP  APGORERT  IS  PIBITE 

C LI°PAPY  ROUTINES  REQUIRED:  NS0 1 A , GAOSSQ,  ODE. 

C L.  N . TarPpTHE  N JANUARY.  1978 

IMPLICIT  F£AL*fl ( A-B , D-H.O- V , X- Y)  , COMPLEX*  16 (C, B,Z) 

C 0 * * D N /SC/  WC.W  (20)  ,BETAP (20)  ,C, Z ( 20) , N , NM, NF 
C C •MON  /CcNSIS/  PI.TWOPI.ZEBO.ZINP.EPS 
PEA L*8  CD  AdS 

C SET  OP  FRCdLEH: 

EPS  ■ 1.D-0 
CALL  INI? 

C COMPUTE  NODES  AMD  WEIGHTS  FOB  PARAMETER  PPOBLEH: 

NPTSQ  ■ 8 

CALL  01  HIT  (RPTSQ) 


C SOLVE  PARAMETER  PROBLEM: 

I PRIM**  * 1 

CALL  SC30LV  (NN.XPRINT) 

C 

C TEST  ACCUPACY  OP  SOLUTION: 
CALL  TEST 


C DP A W CON 

TOUR  PLOT 

OP  5OL0TIOV: 

CALL 

FITCON 

100  CC NT 

TNUE 

STOP 

1 

END 

//GO. SYSIN 

DD  • 

7 

.0 

.0 

2. 

0. 

99. 

2. 

0. 

-.5 

1.  07 0 

1.070 

-1. 

2 

-2. 

-.5 

-.  2 

-1. 

99. 

. 7 

-2.5 

9 1. 

• 0 

-2.7 

99. 

N 

VC 
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pimn  subroutine  •• 




c«  init 

c 


SUBROUTINE  HIT 

INITIALIZES  CONSTANTS  IN  /CONSTS/  AND  PIOBLIN  DEFINITION 

para.ieters  in  /sc/,  data  roi  tee  geceetbt  op  tee  peoples 

IS  1EAD  IN  PROS  DNIT  5. 

IEPLIC IT  REAL*8IA-B.D-N,C-F,X-T),  CCRPLEX* 16 (C. N,Z) 

LOGICAL  FINITE 
CCRPLEX* 16  DCEPLK 

CORSON  /SC/  NC,B(20)  . BETAS  (JO)  ,C,Z(20)  ,N,RH,NP 
COSEON  /CONSTS/  PI .TNOPI, 1EBO. t IBP , (IS 
COEBON  /GEOfl/  XPIX(20)  .BEAT  (20)  .NCOBI 
DATA  S SNARE  /•IBITV 
CALL  ENTER  (SBNABB) 

SE1  CONSTANTS: 

PI  • 3.1*159  26535  89793  23  DO 
TNOPI  * PI  • 2.  DO 
ZEBO  • (0.D0.0.D0) 

ZINF  ■ ( 1. 070, 1.070) 

READ  INPOT  PARARETEBS: 

READ  (5,201)  B 
NR  > N- 1 
NP  • N ♦ 1 

Z (N>  • (1.00,0.00) 

READ  (6,202)  NC 

R3AD  (6,203)  (N  (R)  , BETAB  (R)  ,E-  1 ,B> 

CCEPUTE  ANGLES  AS  REOOIRED  (NRERE  VALUE  INPOT  IS  99.0): 

DO  1)  R * 1 , N 

IF  (0ETA8  (E)  .BE. 99. DO)  GCTO  10 
KR  * ROD (R»N-2,N) *1 
KP  « fl  JD  (E  , N)  ♦ 1 

JtTA.KKl  * D:rAG(CDLOG((N(EE)-N(E))/(N(EP)-H(KI|)|/PI  - 1 . DO 
IF  (6ETAR  (R)  .LS.-1.D0)  BETAE(R)  * BETAB(R)  ♦ 2. DO 
10  CONTINOE 

CHECK  EOF  VARIOUS  INPUT  EPRORS: 

308  • 3. DO 
DO  IN  • 1 , N 

1 SUB  * SUR  ♦ BETAS  (R) 

IF  (DABS (SUR*2.D0) .LT.EPS)  GOTO  2 
NAITJ  (6,3011 
STOP  2 

2 IF  (FINITE  (V  ( 1)  ) ) GOTO  3 
WRIT!  16,  302) 

STOP  2 

3 IF  (FINITE  (N  (N)  ) ) GOTO  A 
UNITE  (6,30  3) 

STOP  2 

* IF  (SETAR  (NR)  .Nr.O.  DO)  GCTO  6 
N9ITS  (6,  30*1 

5 IF  (3ETAH  (NR)  . BE.  1.  DO)  GCTO  20 
WRITE  (6.306| 

STOP  2 

ONTSRRINE  RUBBER  OF  BOUNDARY  CORPONENTS.  ETC.: 

PASS  1:  ONE  FIXED  POINT  FOR  EACN  INFINITE  VERTEX: 

20  NCORP  « 0 

DO  21  R • 2, NR 

IF  (FINITE  (N  (E)  ) ) GOTO  21 
NCORP  • NCCBP  ♦ 1 
EFIX  (NCORP)  ■ E - 1 
IF  (NCORP.  E0.1)  EFII(NCOEf)  • 1 

21  CONTINOE 

IF  (SCOEP.GT.O)  GOTO  22 

NCORP  ■ 1 

EFIX  (NCORP)  • 1 

PASS  2:  ONE  RATIO  POE  EACE  LINE  SEGRENT: 

22  CONTINOE 

NEO  • 2*NC0NP 
DO  23  E * 1 , BE 

IF  INEQ.EO.BB)  GOTO  30 


IP  (• 

NOT. 

FINITE  (N  (t|)  . 

mg  ■ 

NEO 

♦ 1 

ERAT  (NEO) 

• R 

23 

CONTINOE 

30 

CALL  EXIT 

RETORN 

201 

FOREAT 

(13) 

20  2 

FORRAI 

(2FB 

.0) 

20  3 

FORRAT 

(2D8 

.O.PN.O) 

301 

PjRNAT 

(/’ 

• •• 

ERROR 

If 

302 

FORRAT 

(/* 

• •• 

ERROR 

IN 

)33 

FOREAT 

(/• 

• •• 

ERROR 

IN 

30* 

FORRAT 

(/* 

• •• 

NANNING 

30  5 

FOREAT 

(/' 

• •• 

ERROR 

IN 

KD 
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on  r>  n o 


C»*«  •••••«••  

C«  QI HI? 

SOBROUTIBE  QIBIT(BPTS) 

c 

: cohpotes  bodes  bed  heights  rot  gauss-jacobi  quadrature 

c 

IRPLICIT  REAL«8(A-B,D-H,0-f ,S-TI , COUPLER* 16 (C.H, S) 
LOGICAL  PXHXTE 

CDHROR  /SC/  HC,H(20|  .BETAH  (20)  ,C,Z(20)  .H.HH.BP 
COHBOB  /QUAD/  QBODES ( 32,2 1) ,QBTS (32, 3 1) , HPTSO 
DIBEBSIOB  QESCB(2),  QSCB(}2) 

DATA  59BAHB  /'QIBITV 
CALL  EATER  (SBBARE) 

WRITE  (6,201)  BPTS 
C 


PRIRART  SOBROOTIRE  •• 


BPTSQ  * BPTS 

C 

C P0&  EACH  PIBITE  PERTES  B (R) , CORPUTE  BODES  ABD  HEIGHTS  POE 
C ORE-SIDED  GAUSS-JACOBI  QUADRATURE  AIOBG  A CORPS  BEGIHHIIG  AT  Z (K) : 
D1  11  > 1 , B 

1 IP  (PIBITE  (B  (R) ) ) CALL  GAUSSQ (5, HPTSQ.O. DO. BETAH (R) ,0, 

E QESCS.QSCR, QBODES  (1.R)  .QBTS(I.R)) 

C 


0 DEPUTE  NODES  ABD  HEIGHTS  POR  PURE  GAUSSIAB  QUADRATURE: 

CALL  GAUSSQ (5, BPTSQ, 0*D0f0. DO , 0 , QESCR , QSCR , QBODES ( 1 , BP)  , 
E QBTS  ( 1 , HP)  ) 

C 

CALL  EXIT 
S2TUBB 

C 

201  PORRAT  (•  BPTS  -’.IS) 

EBD 


C»  TEST  PRIHART  SOBROUTIBB  •• 


SUPROOTIBE  TEST 

TESTS  TR!  COHPUTED  HAP  POR  ACCOIACT. 

IRPLICIT  REAL*8(A-B.D-H,0-P,Z-T) , COH PLUS* 16 (C, B, Z) 

REAL*9  CDABS 
LOGICAL  PIBITE 

COHROB  /SC/  BC,B(20)  .BETAHI20)  ,C,Z(20),B,RB,BP 
COHROR  /COBSIS/  PI .TBOPI, ZE BO, ZIBP, EPS 
DATA  SBBARE  /‘TEST’/ 

CALL  EBTER  (SBBARE) 

T5ST  LEBGTR  OP  RADII: 

RADE1X  • 0. DO 
DO  10  R * 2,R 

IP  (PIBITE  (H(E)))  BADE  • CDABS  (BC  - BSC  (ZERO,  Z (R)  , B (R)  , R) ) 
IP  (.ROT.  PIBITE  (B(R>)>  BADE* 

6 CDABS  (RSC<(.  IPO,.  IDO)  ,Z  (R-1)  ,H(E-1)  ,R-1> 

6 - HSC(|.  IDO,.  IDO)  ,X(R*1)  ,H(R«1),E»1|) 

RADERS  • DHAX1 ( t ADBH X , RADI) 

10  COBTIBOE 

BRITE  (6,201)  RADERS 
C 

CALL  EXIT 

RETORB 

C 

201  PORRAT  (/•  RADERS:  * , D 12.  B) 

ERD 
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r>  o oooonooo«io 


C»  SCSOLT 


PI  TRAIT  SUBROOTIII  •* 


l - 

SOBIOOTIIE  SCSOLT (IH.IPIIIT) 

C 

C THIS  S'JBROOTIIE  CCHPOTES  TIE  ICCESSOIT  PIIIBETEIS  C AID  Z ( K)  . 

C THE  PB0BLE9  IS  BOLTED  BT  PIIDIIG  THE 

C SOLUTION  TO  A SISTER  OP  1-1  BOIUIEAI  EOUATIOIS  II  TRE  1-1 
C ONXNOHIS  I (11  ...  ..T  (1-1)  , RBICH  AIE  IELATED  TO  TRE  POUTS 
C Z (R)  at  TRE  POIRULA: 

III)  ■ 108  ((TR(R)-TR(R-1))/(TH(R»1)-TR(I)M  (1) 

THEBE  TR(I)  DEIOTIS  THE  AISOBEITOP  Z|R). 

SUBROUTINE  SCPUI  DEPIIBS  THIS  SISTER  OP  ICUATIOIS. 

THE  ORIGINAL  PIOBLEB  IS  SUBJECT  TO  TIE  CCITIAIITS  TRIE)  < TH(N*1), 
3UT  THESE  VANISH  II  TRE  TIAISPOIRATICN  PIOH  Z TO  I. 

SEE  RA1I  PROG  EAR  POP  PORTREI  COHSEITS. 


IMPLICIT  REAL»8(A-B,D-H,0-V,X-I),  COR  HEX*  It  (C,  I . Z) 

C3SH0N  /CCNSTS/  PI , THOPI,  ZE  K , X I IP,  E f S 

OTHERS  101  AJIIT  (20,20)  , SCR  (TOO)  , MALI  19)  , 1(19) 

EXTERNAL  SCPOI 

DATA  SEIARE  /'SCSOLT'/ 

CALL  ENTER  (SBRAPE) 


„ INITIAL  iUZSS  (YTBTICES  EQUALLY  SPACED  APOOID  CIRCLE) i 
DO  ) I ■ 1,IR 
3 T (R)  ■ 0. DO 


C ISO  11  COITIOL  PAIAFETBIS: 

DSTEP  ■ 1.D-8 
DHA1  • 1.D1 
ACC  • EPS 

•Afrai  ■ in  • is 

c 

C S )L  VE  IOILIIEAP  SYSTER  IITH  ISOU: 

CALL  NS01A (IR,Y, Ff  AL, A Jllf # DSTEP ,DH AX , ACC, RAX rON, I PA  I NT, SCI, SCPOI) 
CALL  YZTRAI  (T) 


C PAINT  PESULTS: 

CALL  SCOOTP 


CALL  EXIT 
BETOII 

END 


C*  YZTRAI  SOBCRCTIATE (SCSOLf ) SOBROOTIIB  •• 


SUBROUTINE  YZTRAI (T) 

r 

C lAANSPOMS  Y(K)  TO  Z(K).  SEB  CORREITS  IN  SOBROOTIIB  SCSOLY. 

c 

ISPLIC IT  REAL*8(A-B,D-R,0-?,X-Y) , CORPLEX*  16 (C,l,t) 
CORPLBI* 16  DCRPLX 

COBHOI  /SC/  1C,  I (20)  , BETA R (20)  ,C, Z ( 20) , I , Ifl, IP 
COHHOI  /CCISTS/  PI, T 10 PI, XERO,ZIIP,tFS 
DIHEISIOI  Y ( 1) 

C 

DTK  • 1. DO 
TRSOR  • DTB 
00  1 K • 1,IH 

DTR  « DTB  / DEXP (T  (R) ) 

1 THSOB  • TRSOR  ♦ DTR 
C 

DTR  ■ TIOPI  / TBSOR 
TRSOR  ■ DTB 

Z(1)  • DCRPLX  (DCOS  (DTB)  ,DSIB(DTB) ) 

DO  2 R • 2,IR 

DTR  ■ DTB  / DEXP  ( Y (B- 1)  ) 

TRSOR  • TRSOR  ♦ DTB 

2 Z (K)  ■ DCRPLX (DCOS  (TBSOB) ,DSIN  (TRSOB) ) 


RETURN 

BID 


SO BORDIIATE (SCSOLT)  SOBRODTIIE  •• 


• sc  for 

SOBR  3UTIIE  SCPOI  (IDI H , T,  PTA  L) 

THIS  IS  THE  FOECTIOI  (ROSE  ZEtO  RUST  BE  FOOIO  II  SCSOLT. 

I RFLIC IT  BEAL'S (A-B, D-H,0-T,I-T) , COHPLEI' 16 (C.»,Z| 
REAL'S  CDABS 
LOGICAL  FIBITt 

DIRIISI3I  FFIL(IDIB)  ,T  (1DIH) 

COR  HOI  /SC/  »C,B (ZO) ,BET»H (20) ,C,Z(20) ,1,11,11 
CORIOI  /CCISTS/  PI,TTOPI,ZEB.ZIIP,IIS 
CORIOI  /GEOH/  KPIK20)  .HIT  (20)  .ICOHF 

TRAMS  FOBS  ((E)  TO  Z (I)  : 

CALL  TZTRAI  (T) 

SAT  OF:  CC  SPOTS  IITEGIAL  FlOfl  0 TO  Z (It S 
HOIIOI  • ZOOAD(ZERO,O.Z(I|  . D 
C » (I  (I)  -1C)  / TDETOH 

CASE  t:  1(1)  AID  T(K'I)  FIIITE: 

(CC  SPOTS  IITEGIAL  ALOIG  CHORD  Z (Kl -Z  (K*  1)  I : 

1FIR5T  - 2'ICOHP  ♦ 1 
IF  (HFIRST.GT.IR)  GOTO  11 
DO  13  NEQ  • MPIRST , T H 
KL  - KRAT(IEO) 

KR  • 11*1 

ZIIT  • Z0OAD(Z(KL)  ,KL,Z  (K  D .11) 

FVAL(ISQ)  - CDABS  (I  (Kl) -I  (HI)  ) - CDABS  (C'ZIIT) 

13  C3KTIIUE 

CASS  2:  1(1*1)  IIPIIITE: 

( COMPUTE  CCITODR  IITEGIAL  ALCIG  RADIOS  0-Z(R)|: 

1 1 0 3 23  If E IT  • 1.IC0RP 
KR  • KFIZ (ITER T) 

ZIST  • ZOOAD (ZERO, 0( Z (KR)  «KF) 

ZPTAL  « I (RE)  - *C  - C'ZIIT 
FTAL (2'*TSRT-1)  • DREAl  (ZFTAL) 

FTAL  (2'ITERT)  • DIRAG(ZFTAL) 

20  COITIIOE 
BSTOII 


EIO 

A'OOTP 


3II80RDIIATE  (SCSOLT)  SOBROOTIIE  •• 


3JBR0UTIIE  SCOOTP 

PRUTS  RESULTS  (TARIABLES  II  CORRCI  BLOCK  /SC/) 

ISPLIC IT  REAL'S ( A-B , D-R.O-T ,1- T) , COHPLEI' 16 (C.l, Z) 

LOGICAL  FIIITE 

CORHJI  /SC/  1C, I (20)  .BETAS  (20)  ,C,Z  (201  ,I,IR,IP 
C0RR3I  /COISTS/  PI,TTOPI,ZEBO.ZIIP,IIS 

TRITE  (6,102) 

DO  1 K • 1,1 

THDPI  • DIRAC  (CDLOC(Z(K)t)  / PI 
IF  (TRDPI.LE.O.DO)  TROPI  • THDPI  ♦ 2. DO 

IF  (FI  HITE  (T  (K)  ) ) TRITE  (6,  103)  K , I (K)  , THDPI  , BETA  R (K)  . Z (K) 
1 IF  (.l)T.FIRITEII(K)))  TRITE  (6.10T)  K,  THDPI,  BETAS  (K|  , I (!) 
TRITI  16,106)  TC,C 

RETORI 

102  FORSAT  {//•  RESULTS:'//  

6 • I*. lOI.'SW.UZ.'t*  !»/**•.  i«.'im,c«f. 

6 ISI.'Z(K)'/  . , , 

10|*FOIRAT  IIJ,'  ^ C,F6.  J,  '.FS.J.')  '|F20.  IT, FI*. 5, 

6 31.'  C.FH.  12,',',P1S.12,'>  ')  

10T  FOIRAT  (II,'  IIFIIITT  ' ,F20, IS, FIT, 5, 

r.  31,'  (•.Fli.lR.'.'.FIS.II,')') 

105  FOIRAT  (//•  *C  • (•  ,D22. 15, • ,D22.  15,')  •/ 

A • C • ('  ,D22«  15,  ,022.  15,')  •/) 

RID 
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C*  ZQUAD 


SECONDARY  SUBROUTINE  •• 


FUNCTION  ZQOAD(ZA.KA,ZB,kB) 

C JN  POT! 5 THE  CORFU!  LINE  INTEGRAL  OF  XFNOD  FIOI  Xk  TO  X(  XLONG  k 
STRAIGHT  LINE  SEGBZNT  NITNIN  TNI  UNIT  DISH.  FUNCTION  XQUkDI  IS 
CELLED  TWICE,  ONCI  FOB  EkCB  HkL  F OF  THIS  IBTIGFkL. 

ISFLICIT  REALM  (A-B.D-R.O-T.I-T).  COHFLBX*  It (C.N.X) 

TRIO  • (Xk  ♦ IB)  / 2.  DO 

ZQUkD  > ZQUAD1(ZA,ZHID,kA)  - XOOkD I (XB. XSID, KB) 

B STUBS 
BHD 

• ZQUAD1  SOBOFDIHkTI (ZQUkD)  SUBNOOTIN1  *• 

FUNCTION  ZQOkOIIZk.XB.kk) 

CONFUTES  THE  COHFLEZ  LISE  IRSCUl  OF  ZFBOD  FEOR  Xk  TO  XB  kLOBC  k 
ST) BIGHT  LINE  SECIIRT  NITRIR  THE  UNIT  DISH.  CORFOOND  ONE-SIDED 
GBUSS-JkCOBI  QOBDEBTDBE  IS  USED.  USING  FOBCTION  DIST  TO  DETEEHIRE 
THE  DISTANCE  TO  TRI  REkREST  SINGOLkRITT  X (k)  . 


IHPLICIT  REALM  (k-B.D-R.O-F.X-T)  . COHFLEX*  14  (C.N,  X) 

CCEHON  /COASTS/  FI.TNOFI.XIBO.ZINF.IFS 

REALMS  CDABS 

DATA  IEJFBR  /2.D0/ 

C 

C CHECK  FOE  ZEEO-LENGT!  INTEGRAND: 

IF  (CDABS  (ZA-ZB)  .GT.O.DO)  GOTO  1 

ZQ0AD1  • ZERO 

E3TUER 

<• 

C STEF  1:  ONE-SIDED  GADSS-JACOBI  QUADRATURE  FOE  LEFT  ENDPOINT: 

1 E • DHIN1 (1. DO. DIET (ZA.kk) •HESFRR/CCABS (ZA-ZB) ) 

ZAA  ■ ZA  ♦ P* (XB-ZA) 

ZQUAD 1 • ZQSUH (Sk.Zkk.kk) 

c 

C STEF  2:  ADJOIN  INTERNALS  OF  PURE  GAUSSIAN  QUADRATURE  IF  NECESSARY: 
10  IF  (I.  IQ.  1.  DO)  RETURN 

R • DR  INI  (I.DO.DIST(Xkk.O)  * RES  FIR/C  CAES  (Zkk-ZB) ) 

XBB  • Ekk  • R* (XB-ZA  A) 

ZQUAD1  • ZQUAD  1 ♦ ZQSUH  (ZAA  .ZBE.O) 

ZAA  • XBB 
GOTO  10 
END 

C*««» ••••••• ••••••••••••••••••••••••••••••••••••••«••••••••••••••••• 

C*  DIST  SOIORDIHATE (ZQOA D|  SOBROUTIRB  •• 



c 

FUNCTION  DIST(ZX.kS) 


' 


4 


C DETERSINIS  THE  DISTANCE  FROR  XZ  TO  TRI  NEAREST  SINGULAR! TT  X(k) 

0 OTHER  TRAN  X(RS)  . 

ISFLICIT  REALM (A-B. D-R.O-N ,X-T) , CORPLEX* 16 (C.N, I) 

COPH )N  /SC/  NC.N(20)  , BETA H (20) .C.X(20) .R.NH.NF 
REALM  CDABS 

C 

DIST  • 94. DO 

DO  1 k • 1 , B 

IF  (k.EQ.KS)  GOTO  1 

DIST  • DRIR1  (DIST, CDABS  (ZS-I(k)  ) ) 

1 CONTINUE 
IETURB 
END 

C***« ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• 

C»  IQ JON  SUBORDINATE (SQUAB)  SUBROUTINE  •• 

€••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• 

c 

FUNCTION  EQSUR (SA.SB.RA) 

C 

C CONFUTES  THE  IBTZOEAL  OF  SFROD  FROB  XA  TO  IE  BT  AFFLTINO  A 
C DSX-SIDBO  GOARS- JACOBI  PORROlk  RITE  ROSSIBLB  SIEGOLAEITT  AT  Sk. 

C 

ISFLICIT  BZAlM(A-B,D-B,0-V,X-I),  COEILEt*  1«  (C.R,  E) 

COS  BCR  /SC/  NC,B(20)  . BETA  E (20)  ,C,E  (20)  .R.RN.IF 
CONNOR  /CORSTS/  FI.TROFI, EERO, SIRF, IRS 
CORRDB  /QUAD/  QBODES  1(72)  , QfTS  (612)  .ERTSQ 
MALM  CDABE 
C 

XS  > 1 110 

XR  ■ ( XI- XA)  / 2. DO 
XC  ■ (IAMB)  / 2. DO 
X • RA 

IF  (k. IQ.O)  k ■ I* 

11  • JJMK-1)  ♦ 1 

12  • II  ♦ RRTEQ  - 1 


PRIMARY  SU6R0UT  INfc  •• 


FUNCTION  WSC(Z7tZO,kOtKZO> 

C 

C iNTfcGHATtS  F H 0 H ZO  TO  ZZ  TO  CCMPUTt  N V ALUc  CORN fcS PONDING  TO  ZZ 
C 

IMPLICIT  NtAL*0( A-B,0-H,Q-V ,X-t ) , COMPLEX* 16<C,W ,Z) 

COMMON  /SC/  WC,N(20)  , BETAM< 20 ) # C , 2 < 20 ) ,N,NM,NP 
C 

NSC  « WO  ♦ C • ZLUAD(Z0,KZ0,ZZ,0) 

C 

RETURN 

END 


c*  zsc 




PlinilT  SUBROUTINE  •• 


FUNCTION  ZSC  (NN,Z0,N0,RZ0) 

C 

C IMPUTES  Z(VN).  FIRST  ODE  IS  CELLED  TO  GET  Al  INITIAL  E5TI SATE; 
Z PH  E S ZNEMT  IS  CALLED  TO  GET  THE  PINAL  ARSIER. 

C 

IMPLICIT  REAL*8 (A-B,D-H,0-?,X-l) , COMPLEX*  16 (C, N, Z) 

DIHERSTON  SCR  (M2),  ISCB(S) 

EXTE9NAL  ZPODE 

CCPHJH  /SC/  WC,W (20)  , BETAS (20)  ,C,Z  (20) ,N,NH,HP 
COP HON  /Cf  N STS/  P I , T VOPI , ZE NO , Z I I P, Z F S 
COHHM  /7SCC0P/  CDNDT 

C 

C 35T  INITIAL  GUESS  Z1  VIA  ODE: 

Z1  ■ ZERO 
P ■ 0.  DO 
IPLAG  « -1 
RELE3R  * 0. DO 
ABSE3B  • S.0-3 
CDNDT  * (NN-NC)/C 

CALL  ODE  f ZPODE, 2,Z1,T, 1.D0, R£LEBP,i§SERR, I FLAG, SCR, I SCR) 

IF  (IFIAG.NE.2)  WRITE  (6,201)  I FLAG 


C SEEINE  4NS2XP  >14  ZHENT : 

C4LL  ZNENT (Z1,NN,EPS,KZ0) 

ZSC  * Z1 

c 

231  f OS  R 4T  (V  •••  NONST4ND4B  D 1ETDII  PPOR  ODE  IN  ZSC:  IEL4S  ••,12/) 
IETORN 
(NO 

C*  

C»  ZEODE  SUBOPDIN4TE  (ZSC)  SUBROUTINE  •• 


3 UBR  )OTXN  E ZPODE (T,ZZ, ZDZDT) 

ORP'JTtS  THE  PBNCTION  ZDZDT  HEEDED  BT  ODE  II  ZSC. 

IBPLICIT  RE4l»8(4-B.D-N,0->,X-T).  CORPLEI*  16 |C,H, Z) 
CORSON  /Z SCCCR/  CONDT 


ZDZDT  • CDNDT  / ZPNODIZZ.O) 

IETORN 

END 


S0BCNDIN4TE  (ZSC)  SUBNOOTIRE  •• 


SU8I30TINE  ZNtNT (Z BOOT. ■■.IBS, RIO) 

ISPLEHBNTS  NENTOB'S  NETBOD  TO  SOLTI  Til  EQOITIOB 
2 (ZBOOT)  • Nl  PON  ZPOOT. 

I1PLICIT  BE4l»B(»-B,D-N,0-»,I-T)  , CON  FLEX*  1*  (C , ( , I) 

COP H J N /SC/  HC,N (20) , BET4R ( 20) ,C,Z (20) .N.NN.NB 

DO  1 ITER  ■ 1,»0 
ZBOOT0  • ZBOOT 

IP  (NZO.ZO.O)  ZPBBT  • HI  - NSCIZBOOTO,  (0.00,0.001  ,NC,0| 

If  (RIO. IE. 0)  1PVHT  • »1  - NSC  (ZROOTO  , X (N10)  , H (RIO)  , BIO) 

ZIOOT  • ZIOOTO  ♦ ZPMT/(C»ZFNOD  (ZBOOTO.O)  ) 

IP  (COBBS  (ZPINT1  .LT.  EPS)  BSTOBB 
1 CONTINUE 

TRITE  (6,201) 

IETUBI 

201  E0NH4T  (/•  •••  ENNON  IN  ZNINT:  1C  CCITIICINC1  II  10  ITEIBTIOIS1) 
BID 


SECONDARY  SUBROUTINE  • • 


C*  ZPRJD 

C»* ••••••••• •••••••••••••••••••••••••••• 

c 

poictici  zprod (zz* rs) 

c 

C CJKPUTBS  T HE  INTEGRAND 

c 

C N 

C PICO  (1-ZZ/Z  (K)  ) ••BBTAH(R)  f 

C R«1 

C 

C TARING  AHG01ENT  OILT  (ROT  HODOLOS)  POI  TEBH  K • RS. 

C 

IMPLICIT  FEAL*8(A-B*D-R,0-f *X- Y) * CCRPLBX*  16  (C*N, t) 
NE4L*8  CDABS 

CORN 31  /SC/  MC*R  (20)  ,BBTAH (20)  *C,Z(20) *N*IH,NP 

C 

ZSOR  * (0.00*0.00) 

DO  t I • M 

ZTSP  • (1.00*0.00)  - ZZ/Z  (R) 

IP  (R.EQ.RS)  ZTRP  - ZTHP  /CDABS(ZTHP) 

1 ZSOR  • ZSUH  ♦ BBTAH (R) *CDL0G (ZTRP) 

ZPROD  • CDSIP (ZSOR) 

&ETORR 

END 


C*  PINITB  SECONDARY  SUBROUTINE  •• 

c 

PORC-DN  PINITB  (Z| 

C 

C iZIOFNS  TBUE  IP  ABO  ORLY  IP  Z IS  ROT  IBPIRITE 

C 

IMPLICIT  PEAL* 8 (A-B,D-H*0-?*I-Y)  * CC fl PLEX*  IS  (C * N , Z) 

LOGIC! t PINITB 

TOPRDN  /CONSTS/  PI * TROPI , ZB  BO , III P* IP S 

n 

FINITE  - DFEAL  (Z) .HB.DBEAL  (71 IP) 

RETORN 

BID 


C*  ENTER  SECONDARY  SOBBOOTIBB  •• 


5UBF  30 TIN  B ENTER (SBNARE) 

C 

C STARTS  TIRING  TIRE  SPENT  IB  SOBBOOTINB  NITH  NAH8  SBMARB. 

C 

IflPLICII  BEAL* 8 (A-B,D-R*0~?*X-Y) * COR PLBX*  16  (C , R* Z) 

OEBOB  /TIBI/  TBBTBB 

C 

CALL  LBPT1A  (TENTS!) 

MBITS  (6*201)  SBNARE 
BSTOII 

201  P3BHAT  (//II, 80  (*!•)  * • BITE  BUG  ».A8) 

EBD 

••••••••••••••••••••••••••*••*••••••••••••••••••••••••••••••••* 

c*  EXIT  SECOROABY  SOBBOOTIBB  •• 

CMt«M*MM*t«MM*M*M***t**Mi*MM*«***M***t*M***M*****»**** 

c 

SOBB  30TIBB  HIT 

C 

C PRUTS  TIBI  SPBIT  IB  SOBBOOTIBB. 

C 

IRPLICIT  RIAL*8(A-8,D-R,0-f *I-Y) , CORPLEX*  16 (C*B*t) 

CORROB  /TIRE/  TBBTEB 

C 

CALL  LBPT1A (TEIIT) 

TIRE  - TBBTBB  - TEIIT 
MBIT!  (6,201)  TIRE 

BBTURB 

C 

201  P 3IRAT  (H*  80  (•!'),'  TIBE  IUPSID;  • , P7.  3*  • SECS.*/) 

IRD 


45 
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